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ON  OPTIMAL  LINEAR  SMOOTHING  THEORY 


J.  S.  Meditch* 

ABSTRACT 

The  algorithm  for  generating  the  smoothed  estimate  x(t/t  +  T)  of  the  state 
x(t)  of  a  continuous  linear  system,  where  t  is  continuous  time,  T  is  a  positive 

real  constant,  and  t  +  T  is  the  time  of  the  most  recent  measurement,  is  developed. 
A  linear  matrix  differential  equation  whose  solution  gives  the  covariance  matrix 
of  the  smoothing  error  x(t|t  +  T)  «  x(t)  -  x(t|t  +  T)  is  then  derived.  Computa¬ 
tional  aspects  involved  in  mechanizing  the  algorithm  are  discussed  in  terms  of 
the  algorithm's  dependence  on  the  solution  of  the  prediction,  filtering,  and 
fixed-point  smoothing  problems.  The  results  are  then  discussed  in  terms  of  the 
classical  Wiener  smoothing  problem. 


Department  of  Electrical  Engineering,  Northwestern  University,  Evanston,  Illi¬ 
nois  60201.  This  work  was  supported  in  part  by  the  U.  S.  Army,  Navy,  and 
Air  Force  in  a  Joint  Services  Electronics  Program  under  Office  of  Naval  Research 
Contract  Number  N00014-60-C0020-A03. 
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l. 

1.0  Introduction 

The  continuous  optimal  linear  smoothing  problem  consists  in  determining  a 
minimum  mean-square  error  estimate  x(t|T),  t^  *  t  ^  T ,  of  a  state  vector  x(t) 
based  on  measurements  z(<7)  given  over  the  entire  interval  t  ^  a  *  T  where 

;  -  F(t>  x  +  w(t>  CD 

z(t)  i  H(t)  x(t)  +  v(t)  (2) 

In  Eqs.  (1)  and  (2),  x  is  an  n-vector,  the  state;  z  is  an  m-vector,  the  measure¬ 
ment;  w  is  an  n-vector,  the  disturbance;  v  is  an  m-vector,  the  measurement  error; 

F  and  H  are  continuous  n  x  n  and  m  x  n  matrices,  respectively;  t  denotes  time;  and  the 

dot  denotes  the  time  derivative.  We  assume  that  w(t)  and  v(t)  are  independent 
Gaussian  white  noise  processes  with  identically  zero  means  and  covariance  matrices 

E[v(t)  v'(<T)  ]-  Q(t)  6(t  -  a) 

E[  v(£)  v'(<J)  ]  -  R(t)  6(t  -  o) 

for  all  t  and  cj  where  E  denotes  the  expected  value,  the  prime  denotes  the  trans¬ 
pose,  6  is  the  Dirac  delta  function,  Q(t)  is  a  symmetric  n  x  n  positive  semidefinite 
matrix,  and  R(t)  is  a  symmetric  m  x  m  positive  definite  matrix. 

We  assume  further  that  the  initial  time  t  is  fixed,  and  that  x(t  )  is  a  zero 
mean  Gaussian  random  n-vector  which  is  independent  of  u(t)  and  v(t)  for  all  t  and 
whose  covariance  matrix  P(tQ)  *  E[  x(tQ)  x*(t  )]  is  a  symmetric  n  x  n  positive  semi¬ 
definite  matrix. 

We  define  the  smoothing  error  by  the  relation 
X(t|T)  .  x(t)  -  X(tjT) 
and  the  mean-square  smoothing  error  by 

S  -  E[  x'(t|t)  x(t|r)  ] 

where  t  *  t.  We  call  an  estimate  x(t|T)  that  minimizes  S  an  optimal  smoothed 
estimate. 

We  note  that  there  are  basically  three  separate  cases  to  be  considered  in  the 
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smoothing  problem: 


1.  Fixed-Interval  Smoothing.  In  this  case,  the  interval  [t^.T]  over  which  the 

measurements  are  given  is  fixed,  and  we  obtain  an  optimal  smoothed  estimate 

of  x  for  all  t  e[t  ,t ].  This  case  has  been  treated  by  Bryson  and  Frazier^', 
0  (2) 

and  Rauch,  Tung,  and  Striebelv  .  The  results  are  of  particular  significance 
in  post-experimental  data  analysis  where  one  needs  to  obtain  a  ’’refined" 
estimate  of  the  state  vector  of  a  physical  system  over  the  system's  entire 
operating  time. 


2.  Fixed-Point  Smoothing.  If  an  optimal  smoothed  estimate  of  x  is  desired  at 

only  one  value  of  t  6[t  ,T],  we  call  x(t|T)  a  fixed-point  smoothed  estimate. 

Such  an  estimate  can  be  obtained  using  fixed-interval  smoothing,  but  the 

(3  41 

procedure  is  computationally  inefficient  for  this  purposev  *  .  Moreover, 

if  the  terminal  measurement  time  T  is  not  specified  a  priori,  as  might  be  the 
case  in  an  "on-line"  smoothing  problem,  fixed-interval  smoothing  is  not 
applicable  here.  The  algorithm  for  optimal  fixed-point  smoothing  for  continu¬ 
ous  linear  systems  which  begins  with  the  optimal  filtered  estimate^  of  x  at 

the  fixed  time  t  and  "updates"  this  estimate  recursively  as  more  measurement 

(41 

data  become  available  has  been  developed  by  Meditchv  ' .  The  procedure  is 
applicable  to  both  post-experimental  data  analysis  and  "on-line"  data  pro¬ 
cessing  problems  where  one  requires  a  smoothed  estimate  of  a  physical  system's 
state  at  some  critical  time  during  the  system's  operation. 


3.  Fixed -Lag  Smoothing.  Now  suppose  that  T  is  replaced  by  t  +  T  where  T  *  constant 

>  0  and  t  is  variable  with  t  *  t  .  Then,  we  see  that  x(t!t  +  T)  is  a  "running" 

o 

smoothed  estimate  which  "lags"  behind  the  time  of  the  most  recent  measurement 
by  a  fixed  amount  T.  For  obvious  reasons,  we  call  x(t|t  +  T)  a  fixed-lag 
smoothed  estimate.  Such  estimates  are  primarily  of  interest  in  communication 
and  telemetry  systems  where  one  wishes  an  "on-line"  smoothed  estimate  of  the 
state  of  the  message  or  data  transmitted.  The  intuitive  justification  for  the 
lag  is  that  we  would  expect  less  mean-square  error  in  the  smoothed  estimate 
than  in  the  predicted  or  filtered  estimates  x(t|o),  t  >  o  and  t  ■  cr,  respec¬ 
tively.  * 


In  this  paper,  we  shall  develop  the  algorithm  for  optimal  fixed-lag  smoothing  for 
continuous  linear  systems  of  the  type  described  by  Eqs.  (1)  and  (2).  We  shall 
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also  derive  the  matrix  differential  equation  whose  solution  gives  the  covariance 
matrix  of  the  fixed-lag  smoothing  error 

x(t|t  +  T)  -  x(t)  -  x( t | t  +  T) 

Our  approach  consists  in  considering  the  limiting  case  of  the  tixed-lag 
smoothing  solution  for  discrete  linear  systems  wherein  the  time  between  measure¬ 
ments  is  made  arbitrarily  small.  The  equations  for  discrete  optimal  fixed-lag 

(3) 

smoothing  are  well-known'  ,  and  will  be  used  as  the  starting  point  in  our  work. 
The  limiting  process  that  we  shall  utilize  is  due  to  Kalman^. 

2.0  Optimal  Fixed-Lag  Discrete  Linear  Smoothing 

We  begin  by  considering  the  discrete  system  analog  of  Eqs.  (1)  and  (2)  which 
can  be  expressed  in  the  form 

x(k  +1)  -  *(k  +  l,k)  x(k)  +  w(k)  (3] 

z(k  +  1)  -  H(k  +  1)  x(k  +  1)  +  v(k  +  1)  (4] 


where  v,  w,  x,  z,  and  H  are  of  the  same  dimensions  as  they  were  in  Eqs.  (1)  and 
(2);  $  is  an  nxn  matrix,  the  state  transition  matrix;  and  k  ■  0,1,...,  is  the 
discrete  time  index.  We  assume  that  w  and  v  are  independent  Gaussian  white 
sequences  with  identically  zero  means  and  covariance  matrices 

E  [w(J)  w'(k)]  -  Q(k)  6 


E  [  v(J)  v'(k)]  -  R(k)  6 


for  all  j  and  k  where  6  is  the  Kronecker  delta,  Q(k)  is  a  symmetrix  nxn  positive 
semidefinite  matrix,  and  R(k)  is  a  symmetric  m  x  m  positive  definite  matrix. 

We  further  assume  that  the  initial  state  x(0)  is  a  zero  mean  Gaussian  random 
n-vector  which  is  independent  of  w(k)  and  v(k)  for  all  k  and  whose  covariance 
matrix  P(0)  -  E  [x(0)  x'(0)]  is  n  x  n  symmetric,  and  positive  semidefinite. 

We  let  N  be  some  fixed  positive  integer  and  denote  a  fixed-lag  smoothed  esti¬ 
mate  of  x(k) ,  given  measurements  up  to  and  including  the  one  at  time  k  +  N,  by 
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x(k|k  +  N)  where  k  *  0,1,...  .  We  define  the  corresponding  smoothing  error  and 
its  mean-square  value  by  the  relations 

x(k|k  +  N)  *  x(k)  -  x(k|k  +  N) 


S  -  E[x'(k|k  +  N)  x(k|k  +  N)] 

respectively.  We  call  an  estimate  x(k |k  +  N)  that  minimizes  S  the  optimal  fixed- 
lag  smoothed  estimate.  It  has  been  shown'  '  that  this  estimate  is  given  recur¬ 
sively  by  the  system  of  n  first-order  difference  equations 


x(k  +  1 |k  +  1  +  N)  =  $(k  +  l,k)  x(k|k  +  N) 


+  C(k  +  N,k  +  1)  K(k  +  1  4fl) 


.[z(k  +  1  +  N)  -  H(k  +  1  +  N)  x(k  +  1  +  K|k  +  N)] 
+  Q(k)  i ' (k,k  +  1)  P*1(k) 

.  t  i(k|k  +  H)  -  S(k)]  (5) 


for  k  =  0,1,...,  where 


C(k  +  N,k  +  1) 


i«k+l 


J(i)  -  P(i)  $'(i  +  l,i)  M  (i  +  1) 


K(k  +  1  +  N)  «=  P(k  +  1  +  N)  H*(k  +  1  +  N)  R  (k  +  1  +  N) 


[  ]  denotes  the  matrix  inverse,  and  x(k)  and  x(k  +  1  +  N|k  +  N)  are  the 
filtered  and  predicted  estimates  of  x(k)  and  x(k  +  1  +  N) ,  respectively.  In 
addition,  the  n  x  n  matrices  P  and  M  are  the  covariance  matrices  of  the  filtering 
and  prediction  errors 


»«»  - 
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X(j  +  1)  -  X(j  +  1)  -  X(j  +  1) 


and 

A 


X(J  +  1|J)  *  x(j  +  1)  -  x(j  +  l! j) 

respectively. 

The  filtered  and  predicted  estimates  along  with  their  corresponding  error 
covariance  matrices  are  governed  by  the  set  of  re'iations^*^ 


x(j  +  l|j)  -  *(j  +  1,J)  x(j)  (9) 

x(j  +  1)  -  x(J  +  1|J)  +  K(J  +  1)  [z(j  t  1)  -  H(J  +  1)  x(j  +  1 | j) 3  (10) 

K(J  +  1)  -  M(j  +  1)  H*(j  +  1)  [  H( j  +  1)  M(J  +  1)  H’(j  +  1)  +  R( j  +  l)]*1  (11) 

M(j  +  1)  -  t(j  +  l.J)  P(j)  «'(j  +  l,j)  +  Q(j)  (12) 

P(j  +  1)  -.[I  -  K(j  +  1)  H( j  +  1)]  M( j  +  1)  (13) 


for  j  -  0,1,...,  where  x(0)  -  0,  P(0)  =  E[  x(0)  x'(0)  ],  and  I  is  the  n  x  n  identity 
matrix. 

The  n  x  m  matrix  K  which  is  given  by  either  Eq.(8)  or  Eq.(ll)  and  which  also 
appears  in  Eq.  (5)  is  called  the  optimal  filter  gain. 

In  order  to  initiate  fixed-lag  smoothing,  we  note  that  at  k  ■  0,  x(N  +  1 | N) , 
x(0)  ,  and  x(0|N)  must  be  input  to  Eq.  (5).  The  first  two  of  these  follow  directly 
from  the  results  for  prediction  and  filtering.  However,  x(0|N)  must  be  obtained 
from  the  fixed-point  smoothing  filter  by  starting  with  x(0)  s  0  and  processing  the 
measurements  at  k  *  1,...,N.  The  equations  necessary  to  do  this  are  given  else- 
where v  *  '  for  the  discrete  case.  They  are  of  no  consequence  here  since  we  are 
concerned  with  the  continuous  case.  After  wc  have  developed  the  algorithm  for  con¬ 
tinuous  fixed-lag  smoothing,  we  shall  show  in  detail  what  procedure  must  be  followed 
to  obtain  the  appropriate  initial  conditions.  In  any  event,  it  is  clear  that 
iixed-lag  smoothing  depends  upon  inputs  from  both  the  predictor-filter  and  the 
fixed-point  smoothing  filter. 

The  covariance  matrix  of  the  fixed-lag  discrete  optimal  linear  smoothing  error 

-f 

which  was  defined  earlier  is  given  by  the  first-order  n  x  n  matrix  difference  equation 
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P(k+l| k+14N)  =  M(k+1)  -  C(k+N,k+1)  K(k+1+N)  H(k+1+N)  M(k+l-HI)  C'(k-ffl,k+l) 

-  J_1(k)[  P(k)  -  P(k|k+N)]  j"i'(k)  (14) 


for  k  «  0,1, ,  where 

P(k|k  +  N)  -  E  [x(k|k  +  N)  x(k|k  +  N) ] 

and  all  of  the  other  terms  were  defined  previously. 

The  initial  condition  on  Eq.  (14)  is  P(0|n)  which  must  be  obtained  from  the 
covariance  equation  for  fixed-point  smoothing^*^ .  We  shall  present  the  equation 
whose  solution  will  give  the  initial  condition  for  the  continuous  version  of  Eq. 

(14)  in  Section  4.0 

3.0  Optimal  Fixed-Lag  Continuous  Linear  Smoo thing 

Let  us  assume  that  the  system  of  Eqs.  (3)  and  (4)  has  been  developed  by  dis¬ 
cretizing  the  system  of  Eqs.  (1)  and  (2).  The  corresponding  fixed-lag  smoothing 
filter  is  then  defined  by  Eqs.  (5)  through  (8).  We  now  consider  the  limiting  be¬ 
havior  of  this  latter  set  of  equations  as  the  time  between  measurements  is  made 
arbitrarily  small. 

We  let  the  discrete  time  instants  k  and  k  +  1  be  denoted  by  t  and  t  +A  t, 
respectively,  where  At  >  0.  We  let  the  time  interval  N  be  denoted  by  T  *  con¬ 
stant  >  0,  and,  as  a  result,  see  that  k  +  II  and  k  +  1  +  N  become  t+T  and 
t  +  A t  +  T,  respectively.  In  considering  the  limiting  case,  the  covariance  matrices 
Q(k)  and  R(k)  must  be  replaced  by  Q(t)/At  and  R(t)/At,  respectively.  This  is 
nece.sary  in  order  to  obtain  a  physically  meaningful  description  of  the  Gaussian 
white  noise  as  the  limit  of  the  Gaussian  white  sequence.  In  addition,  a  factor  of 
At  arises  in  the  disturbance  term  so  that  w(k)  is  replaced  by  w(t)At.  As  a  result, 
Q(k)  i^  replaced  by  Q(t)At  in  all  covariance  relations  involving  Q(k) .  The  de¬ 
tails  of  the  justification  for  this  procedure  are  given  elsewhere^*^*^ ,  and  will 
not  be  repeated  here. 

Making  these  substitutions  into  Eqs.  (5)  and  (8),  we  have 

x(t+At| t+At+T)  =  $(t+At,t)  x(t| ;+T)  +  C(t t+At‘,  K(t+At+T) 

.[  z( t+At+T)  -  H( t+At+T)  x( t+At+T |  t+T)  ] 

+  Q(t)  $'(t,t+At)  P-1(t)  [  x(t]  t+T)  -  x(t)]£t 


(15) 


7. 


and 


r.( t+At+T)  »  P( t+At+T)  H'  (t+At+T)  R*1(t+At+T)  At 


(16) 


respectively. 

Since  the  state  transition  matrix  satisfies  the  relations 


*(t,T)  -  F(t)  $(t,T) 


and 


*(t,t)  **  i  for  all  t, 


we  see  that  $(t  +  At|t)  and  +At)  can  be  expanded  in  the  Taylor  series* 

$(t  +  A  t  ,t)  *  I  +  F(t)  At  +  0(At2) 


(17) 


and 


**(t,t  +  A  t)  =  I  -  F*  (t)  At  +  0(At  ) 


(18) 


2  2 
respectively,  where  0(At  )  denotes  terms  of  order  (At)  . 

Substituting  Eqs.  (16),  (17),  and  (18)  into  Eq.  (15),  and  rearranging  terms, 

we  obtain 


x(t+At |  t+At+T)  -  x( 1 1  t+T)  =  F(t)  x(t  |  t+T)  At  +  C(t+T,t+At)  P(t+At+T) 

•  H(t+At+T)  R-1(  t+At+T)  [  z(  t+At+T)  -  H(  t+At+T)  x(  t+At+T  I  t+T)  ]  At 
+  Q(t)  P_1(t)  [x(t|  t+T)  -  5(C)]  At  +  0(At2) 


Dividing  through  by  At  and  taking  lira  At  -*  0,  we  then  have  that 

x(t 1 1  Ah  T)  =  F(t)  x( 1 1 1  +  T)  +  lim  i  C(t  +  T,t  +  At)  P(t  +  At  +  T) 

At  -»  0  1 

•  H' (t+At+T)  R_1(t+At+T)[  z(  t+At+T)  -  H(  t+At+T)  x(  t+At+T  I  t+T)  ]} 

+  Q(t)  P*1(t)[  x(t|t  +  T)  -  5(t)]  (19) 


where  it  remains  for  us  to  evaluate  the  second  term  on  the  right-hand  side. 
From  Eq.  (6),  it  is  seen  that 


8. 


C(t  +  T,t)  =  J(t)  C(t  +  T,t  +£t) 


which  can  also  be  written 


C(t  +  T,t  +  A  t)  -  J  (t)  C(t  +  T,t) 


From  Eq.  (7) , 


J*1(t)  =  M(t  +  A t)  +AO  P*1(t) 


Making  the  appropriate  substitutions  into  Eq.  (12),  we  obtain  the  result 
M(t  +  A  t)  =  $(t  +  A  t,t)  P(t)  *'(t  +At,t)  +  Q(t)  At 

=  [I  +  F(t)At  +  0(At2)]  P(t)  [I  +  F(t)At  +  0(At2)]'  +  Q(t)At 
■*  P(t)  +  [F(t)  P(t)  +  P(t)  F'(t)  +  Q(t)  ]  At  +  0(At2) 

=  P(t)  +  0(At)  (2; 

where  0(At)  denotes  terms  of  order  At. 

Substituting  Eqs.  (22)  and  (18)  into  Eq.  (21)  and  simplifying  the  result,  we 
see  that 


J'l(t)  =[  P(t)  +  0(A t)]  [  X  -  F’(t)At  +  0(At2)  ]  P‘l(t) 


=  I  +  0(At) 


Hence,  Eq.  (20)  can  be  written 


C(t  +  T, t  +  At)  *  C(t  +  T, t)  +  0(At) 


from  which  it  immediately  follows  that 


lim  C(t  +  T, t  +  At)  *  C(t  +  T, t) 
At  -♦  0 


From  the  nature  of  the  matrices  P,  H,  and  k  \  it  also  follows  that 


■***«.' 


? 

* 
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lim  P(t+At+T)  H'(t+At+T)  R_  1  ( t+At+T)  =  P(t+T)  H'(t+T)  R'^t+T) 

At  *♦  0 

From  Eqs.  (9)  and  (7),  we  see  that 

x(t  +  At  +  T|t  +  T)  *  $(t  +  At  +  T,  t  +  T)  x(t  +  T) 

-  [I  +  F(t)At  +  0(At2)]  x(t  +  T) 

*  x(t  +  T)  +  0(At) 

Hence, 

lim  [  z(t+At+T)  -  H( t+At+T)  x(  t+At+T  |  t+T)  ]  =  z(t+T)  -  H(t+T)  x(t+T)  (25) 

At  -*  0 

As  a  consequence  of  Eqs . (23 , (24) , and  (25) , the  second  term  on  the  right-hand 
side  of  Eq.  (19)  can  be  evaluated  as  the  product  of  the  three  limits  in  these  equa¬ 
tions,  and  we  are  led  to  the  result 

x(t | t+T)  =  F(t)  x(t | t+T)  +  C(t+T,t)  P(t+T)  H' (t+T)  R_1(t+T) 

*[>(t+T)  -  H(t+T)  x(t+T)]  +  Q(t)  P_1(t)  [x(t  |  t+T)  -  x(t)]  (26) 

which  specifies  the  optimal  fixed-lag  continuous  linear  smoothing  filter. 

Since 

K(t  +  T)  =  P(t  +  T)  H'(t  +  T)  R*1(t  +  T) 

is  the  n  x  m  gain  matrix  of  the  optimal  continuous  linear  filter^’*’^  ,we  can  also 
write  Eq.  (26)  as 

x(t|  t+T)  =  F(t)  x(t|  t+T)  +  C( t+T, t)  K(t+T)  [z(t+T)  -  H( t+T)  x(t-fT)] 

+  Q<t)  P*1<t)  [x(t|t+T)  -  x<t)]  (27) 


(24) 


where  t  ^  t  . 

o 

We  note  that  mechanization  of  Eq.  (27)  requires  the  following  input  data; 
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A.  The  system  matrix  F(t)  and  the  system  disturbance  covariance  matrix  Q(t). 

B.  1'ne  "gain  times  residual"  term  K(t  +  T)  [s(t  +  T)  -  H(t  +  T)  x(t  +  T)]. 

C.  The  time  history  of  the  n  x  n  smoothing  filter  gain  matrix  C(t  +  T,t). 

D.  The  optimal  filtered  estimate  &(t)  and  the  inverse  of  its  error  covariance 
matrix,  i .e. ,  V  *(t) 

E.  The  initial  condition  x(t  ! t  +  T) 

o 1  o 

Item  A  is  generally  given  in  the  problem  sDCcif ication.  Items  B  and  D  follow 
directly  from  the  optimal  filter  results ^’^which  we  repeat  here  for  convenience: 

x  -  F(o)  x  +  K (a)  [  z (a)  -  H(o)  x  ]  (28) 

K  (a)  -  P(a)  H'(a>  R_1(a)  (29) 

P  ■=  F (a)  P  +  PF'(a)  -  PH'(o)  R*X(o)  H(a)  P  +  Q (a)  (30) 

where  a  ;>  t  ,  x(t  )=0,P(t  )  =  E  fx(t  )  x’(t  )*],  and  the  dot  denotes  the  derivative 
o  o  o  u  o  o  J 

with  respect  to  o •  In  this  connection,  we  observe  that  the  filter  of  Eq.  (28) 

must  process  the  measurements  in  the  interval  [t0>t0+  T  ]  before  smoothing  is 

initiated  since  the  term  K(t  +  T)1  z(t  +  T)  -  H(t  +  T)  x(t  +  T)]is  required 

o  o  o  o  J 

in  Eq.  (27)  to  begin  smoothing.  As  a  result,  we  have  two  time  scales  here: 
the  filter  time  scale  c  and  the  smoothing  filter  time  scale  t  where  t;  =  a  -  T, 
i.e.,  the  smoothing  filter  must,  of  necessity,  "lag'*  the  filter  of  Eq.  (28)  by 
T  units  of  time  in  executing  its  data  processing. 

Because  of  the  relationship  between  the  two  time  scales,  we  see  that  P  * (t) 

A 

and  x(t)  which  are  required  in  Eq .  (27)  can  be  obtained  directly  from  P(a)  and 
xp)  which  must  be  computed  first  in  order  tq  mechanize  Eq  (28)  by  introducing 
a  time  delay  of  T  units  into  the  latter  two  quantities. 

We  summarize  our  discussion  here  in  the  form  of  the  block  diagram  shown 
below  in  Fig.  1. 


Q(t)P  -(t) 


(to!to+T) 


C(t+T,t) 


Delay 


L I 


*<'0> 


Smoothing  filter  block  diagram. 

We  turn  now  to  consideration  of  items  C  and  E.  In  particular,  the  initial 
condition  x(tQJto  +  T)  can  be  obtained  by  utilizing  optimal  fixed-point  continu 

ous  linear  smoothing^  to  process  the  measurements  in  the  interval  [t  ,t  +  T]. 

(41  0  0 

The  required  algorithm  isv  ' 


*(t  |°)  ■  B(cr,t  )  K  (a)  [z(a)  -  H(a)  x(a)  ] 


utMtaft**  on^iakJt va 
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I 


where  a  £  t  .  x(t  It  )  =  0,  and  the  n  x  n  gain  matrix  satisfies  the  matrix  differen- 
o  v  o  o' 

tial  equation 


B(a,t  )  -  -  B(o,t  )[  F(°)  +  Q(°)  P_1(a)  ] 


(32) 


where  B(to>to)  «  I. 

The  term  K(a)  [  z (a)  -  H(0)  x(a)  ]  is  obtained  directly  from  the  optimal  filter 
defined  by  Eqs.  (28)  through  (30),  and  the  fixed-point  smoothing  filter  can  be 
operated  on  the  same  time  scale,  i.e.,  simultaneously,  with  the  filter  of  Eq.  (28) v  . 

At  cj  *  tQ  +  T,  fixed-point  smoothing  is  terminated,  and  the  output  x(tQ| t^  +  T) 
of  this  smoothing  procedure  is  input  as  the  initial  condition  for  fixed-lag  smoothing 
(see  Eq.  (27)  and  Fig.  1).  We  note,  of  course,  that  the  fixed-lag  smoothing  filter 
has  been  inoperative  during  the  interval  tQ  ^  cr  £  t  +  T  while  "waiting”  for  the 
required  initial  condition. 

We  conclude  this  section  by  developing  the  matrix  differential  equation  whose 
solution  is  the  fixed-lag  smoothing  filter  gain  C(t  +  T,t).  From  Eq.  (6),  we  see 
that 


k+l4tf 


C(k  +  1  +  N,k  +  1)  »  I  |  J(i) 

i*k+l 


k+N 


IT 


i»k 


J(i) 


J(k  +  i.  +  N) 


J*l(k)  C(k  +  N,k)  J(k  +  1  +  N) 


which  we  choose  to  express  in  the  form 


C(k  +  1  +  N,k  +  1)  J*X(k  +  1  +  N)  -  J'V)  C(k  +  N,k) 


(33) 


Replacing  kby  t,  k  +  lbyt+At,  k  +  Nbyt+T,  and  k  +  l+  Nby  t+At+T 
in  Eq.  (33),  we  obtain 
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C(t  +  At  +  T,t  +  At)  j'X(t  +  At  +  T)  «  J*l(t)  C(t  +  T.t)  (34) 

From  Eq.  (21),  we  recall  that 

J*l(t)  -  M(t  +  At)  i’(t,t  +  At)  P*l(t)  (21) 

From  the  first  line  in  Eq.  (22),  and  Eqs.  (17)  and  (18),  we  see  that 

! 

* 

\ 

M(t  +  At)  i’(t,t  +  At)  «  4(t  +  At,t)  P(t)  +  Q(t)  $'(t,t  +  At)  At 

2  i 

-  [I  +  F(t) At  +  0(At  )]  P(t) 

+  Q(t)  [I  -  F'(t)At  +  0(At2)]At 

-  P(t)  +[F(t)p(t)+Q(t)]At  +  0(At2) 

Postmult iplying  this  result  by  P~*(t),  we  obtain  the  result 

J‘l(t)  “  I  +  [  F(t)  +  Q(t)  P*l(t)]At  +  0(At2)  (35) 

In  an  identical  manner,  we  also  have 

J*1(t+At+T)  -  I  +[F(t+At+T)  +  Q(t+At+T)  P*1(t+At+T)]  At  +  0(At2)  (36) 

Substituting  Eqs.  (35)  and  (36)  into  Eq.  (34)  and  rearranging  the  result 
gives  us 

C(t+At+T,t+At)  -  C(t+T,t)  -[  F(t)+Q(t)P*1(t)  ]C(t+T,t)At 

-  C(t+At+T,t+At)  [  F(t+At+T)-H)(t+At+T)p‘X(t+At+T)]  At 
+  0(At2) 

Dividing  through  by  At  and  taking  lim  At  0,  we  then  have 

C(t+T,t)  -[  F(t)-K)(t)P*l(t)]  C(t+T, t)  -  C(t+T,t)  [  F(t+T)-K)(t+T)P*1(t+T)]  (37) 

which  is  the  desired  result.  The  initial  condition  for  Eq.  (37)  is  C(tQ  +  T,tQ)  * 

A 

i 
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B(to  +  T.t^}  which  is  simply  the  solution  of  Eq.  (32)  evaluated  at  o  *  t^  +  T.  Tn 
this  case,  we  obse'rve  that  computation  of  C(t  +  T,t)  where  t  *  t^  cannot  be  initiat¬ 
ed  until  both  Eqs.  (30)  and  (37)  have  been  solved  over  the  interval  tQ  ^  o  £  t^  +  T 
to  obtain  P(tQ  +  T)  and  B(tQ  +  T,tQ),  respectively.  We  note,  of  course,  that  P(o) 
is  also  required  for  o  >  tQ  +  T.  Finally,  although  the  gain  matrix  C(t  +  T,t)  can 
be  computed  a  priori,  i.e.,  before  fixed-lag  smoothing  is  initiated,  by  solving  Eqs. 
(30),  (32),  and  (37),  the  fact  remains  that  the  smoothing  filter  of  Eq.  (27)  cannot 
begin  functioning  until  x(tQ!to  +  T)  and  the  term  K(t^  +  T)  [z(tQ  +  T)  -  H(to  +  T) 
•x(tQ  +  T)  ]  are  determined  by  fixed-point  smoothing  and  filtering,  respectively, 

over  the  interval  t  ^  a  £  t  4-  T. 

o  o 

4.0  Optimal  Fixed-Lag  Continuous  Linear  Smoothing  Error  Covariance 

Replacing  k  by  t,  k  +  1  by  t  +  At,  k  +  N  by  t  +  T,  and  k+l+Nby  t  +  At+T 
in  Eq.  (14),  and  substituting  into  this  result  from  Eq.  (16),  we  have 

P(t+Atj t+At+T)  =  M( t+At)  -  C ( t+T , t+At) P ( t+At+T) H ' ( t+At+T)  if 1 ( t+At+T) 

H( t+At+T)  M(  t+ At+T) C 1  ( t+T ,  t+At)  At 

-  J*1(t)  [P(t)  -  P ( t | t+T) ]  J*l'(t)  (38) 

For  the  present,  we  focus  our  attention  on  the  first  and  third  terms  of  Eq. 

(38). 

We  recall  from  Eq.  (22)  that 

M(t  +  At)  -  P(t)  +  [F(t)  P(t)  +  P(t)  F'(t)  +  Q(t)  ]At  +  0(At2)  (22) 

Utilizing  the  expression  for  J~*(t)  as  given  in  Eq.  (35),  we  observe  that 

J*1(t)  P(t)  J*r(t)  »  {i  +  [F(t)+Q(t)P*X(t)  ]At+0(  At2)  }  P(t) 

•{i  +  [F'(t)+P*1(t)Q(t)]  At  +  0(At2)} 

-  {  P(t)  +  [F(t)P(t)+0(t)]  At  +  0(At2)} 

*  f  X  +  [F'(t)+P*1(t)Q(t)]  At  +  0(At2)  } 

-  P(t)  +  [F(t)P(t)-H)(t)  ]  At  +  [P(t)F'  (t)+Q(t)  ]  At  +  0(  At2) 
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-  P(t)  +[F(t)P(t)+P(t)F'(t)+2Q(t)]  At  +  0 (it2)  (39) 

Similarly,  we  have  that 

j"1(t)P(t|t+T)J*l'(t)  -[  P(t|t+T)+[F(t)-K)(t)p'1(t)]  P(t|t+T)At  +  0(At2)} 

*{l  +  [F'(t)+P*l(t)Q(t)]  At  +  0(At2) j 
-  P(t |t+T)+  t F(t)-M)(t)p"l(t)]  P(t|t+T)  At 

+  P(t|t+T)  [F(t)+Q(t)P*1(t)]  '  At  +  0(At2)  (40) 

Combining  Eqs.  (22),  (39),  and  (40),  we  see  that 

M(t+  t)-J*1(t)[P(t)-P(t|t+T)]j‘1'(t)  *  P(t|t+T)+[F(t)+0(t)P*1(t)]  P(t |t+T) At 

+  P(t|t+T)  [F(t)+Q(t)P*1(t)],At-Q(t)At 
+  0(At2)  (41) 

Substituting  Eq.  (41)  into  Eq.  (38)  and  rearranging  terms,  we  have 

P(t+At  |t+At+T)  -P(t  |t+T)  «  [F(t)-K)(t)P*1(t)]P(t|t+T)  At+P(t|t+T)[F(t)-K)(t)P*1(t)]  At 
-  C(t+T,t+At)P(t+At+T)H’(t+At+T)R'1(t+At+T)H(t+At+T) 
•M(t+At+T)C'(t*-T,t+At)  At  -  Q(t)  At  +  0(At2) 

Dividing  through  by  At,  taking  lim  At  “*  0,  utilizing  the  results  in  Eqs.  (23)  and 
(24) ,  and  noting  from  Eq.  (22)  that 

lim  M(t  +  At  +  T)  «  P(t  +  T) 

At  -♦  0 

we  obtain 

P(t  |t+T)  ■=  [  F(t)+Q(t)P*1(t)]  P(t  |t+T)  +  P(t|t+T)  [F(t)+Q(t)P*1(t)]  ' 

-  C(t+T,t)P(t+T)H'(t+T)R"1(t+T)H(t+T)P(t+T)C'(t+T,t)  -  Q(t)  (42) 
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which  is  the  result  sought. 

The  n  x  n  matrix  C(t  +  T,t)  in  Eq.(42)  is  the  solution  of  Eq.(37).  The 

covariance  matrices  P  *(t)  and  P(t  +  T)  are  obtained  from  the  solution  of  the 

filter  error  covariance  relation,  Eq.  (30),  as  described  previously.  Finally, 

the  initial  condition  for  Eq.  (42)  is  P(t  |t  +  T)  which  can  be  obtained  by  solv- 

0  0  (4) 

ing  the  fixed-point  smoothing  error  covar*  ance  equation v  7 

P(to|<J)  =  -  B(a,to)P(a)H,(a)R'1(a)H(cr)p(o)Bl(o,to)  (43) 

over  the  interval  tQ  ^  a  ^  tQ  +  T  where  B(^,tQ)  is  the  solution  of  Eq.  (32)  and 
the  initial  condition  for  Eq.  (43)  is  P(t^), 

Solution  of  Eq.  (42)  then  gives  the  covariance  matrix  of  the  optimal  fixed-lag 
continuous  linear  smoothing  error  X(t|t  +  T)  =  x(t)  -  ft(t | t  +  T)and  It  follows  that 

E  t  x'(t  )t  +  T)  x(t  (t  +  T)  ]  “  trace  P(t  It  +  T) 

Finally,  by  noting  the  definition  of  the  optimal  filter  gain  matrix  K,  we  see 
that  Eq.  (42)  can  also  be  written  as 

P(t  |t+T)  =[  F(t)+Q(t)P*X(t)  ]  P(t|t+T)  +  P(t|t+T)  [F(t)+Q(t)P*1(t)]' 

-  C(t+T,t)  K(t+T)  H(t+T)  P(t+T)  C 1 ( t+T , t)  -  Q(t)  (43) 


5.0  Discussion  of  Results 

For  convenience  of  reference  in  the  discussion  to  follow,  let  us  summarize  the 
results  for  optimal  fixed-lag  continuous  linear  smoothing.  The  smoothing  filter 
equation  is 

x(t|t+T)  =  F(t)  x( t | t+T)  +  C(t+T,t)  K(t+T)  [z(t+T)  -  H( t+T)  x(t+T)] 

+  Q(t)  p'*(t)  [x(t/t+T)  -  x(t)]  (27) 


where 

K(t+T)  «=  P(t+T)  H '  ( t+T)  R*1(t+T) 

and  C(t,t+T)  is  the  solution  of  the  n  x  n  matrix  differential  equation 


•t 
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C(t+T,t)  =  F(t)+Q(t)P*1(t)  C(t+T, t)  -  C(t+T,t)  F(t+T)H)(t+T)P*1(t+T)  (37) 

In  these  three  equations  t  ^  t^,  and  the  initial  conditions  x(tQ|  tQ  +  T)  and 

C(t©  +  T,t^)  are  obtained  from  the  solution  of  the  optimal  fixed-lag  continuous 

linear  smoothing  problem  over  the  interval  [  t  ,t  +  T"|.  The  n  x  m  matrix  K(t  +  T)  is 

u  o  o  ** 

the  gain  matrix  for  optimal  continuous  linear  filtering.  The  matrix  C(t  +  T,t)  is 
termed  the  optimal  smoothing  filter  gain. 

The  fixed-lag  smoothing  error  covariance  matrix  equation  is 

P(t|  t+T)  =[F(t)+Q(t)p”^(t)  ]  P(t|  t+T)  +P(tjt+T)  [F(t)+Q(t)P_1(t)  ]' 

-  C(t+T,t)K(t+l)H(t+T)P(t+T)C'(t+T,t)  -  Q(t)  (43) 


for  t  ^  t^  where  the  initial  condition  P(^o|tQ  +  T)  is  obtain  from  the  optimal 
fixed-lag  smoothing  solution. 

Perhaps  the  most  striking  feature  of  the  fixed-lag  smoothing  filter  described 
by  Eq.  (27)  is  that  it  contains  two  "correction"  terms  in  addition  to  the  "homo¬ 
geneous"  term  F(t)  x(t|t  +  T) .  This  is  in  contrast  to  the  familiar  Kalman-Bucy 
filter  of  Eq.  (28)  which  possesses  a  single  "correction"  term  in  addition  to  the 
homogeneous  term. 

The  first  "correction"  term  in  Eq.  (27)  is  a  weighting  of  the  "gain  times 
residual  term"  K(t  +  T)  [  z(t  +  T)  -  H(t  +  T)  x(t  +  T)]  found  in  the  Kalman-Bucy 
filter.  The  function  of  the  smoothing  filter  gain  C(t  +  T,t)  is  to  weight  the 
"information"  in  K(t  +  T)  £z(t  +  T)  -  H(t  +  T)  x(t  +  T)]  and  "reflect"  it  into 
x(t  +  T).We  recall  here  that  the  estimate  lags  the  measurement  by  T  units  of  time. 

The  second  "correction"  term,  on  the  other  hand,  involves  a  weighting  of  the 
difference  between  the  fixed-lag  smoothed  estimate  and  the  filtered  estimate,  both 
at  time  t.  We  can  view  this  difference  as  a  smoothing  vs.  filtering  "error"  signal. 

Now  let  us  recall  that  Q(t)  is  the  covariance  matrix  of  the  system  disturbance 
w(t)  in  Eq.  (1)  and  that  P(t)  is  the  covariance  matrix  of  the  filtering  error  x(t). 

We  note  immediately  that  if  there  is  no  system  disturbance,  then  Q(t)  =  0  for  all 
t  ^  tQ  and  the  second  "correction"  term  vanishes.  This  is  plausible  for  the  follow¬ 
ing  reason.  If  Q(t)  =  0  for  all  t  ^  tQ,  the  uncertainty  in  x(t)  is  due  entirely  to 
the  uncertainty  in  x(tQ),  the  initial  state.  This  uncertainty  can  only  be  reduced 
by  examining  the  measurements  z(t)  for  t  ^  tQ  in  which  case  the  difference  x(t|t  +  T) 
-  x(t)  contains  no  "new  information"  not  already  present  in  the  residual  z  -  H  x. 

Now  suppose  Q(t)  ^  0  for  t  ^  tQ.  Then, if  x(t)  is  an  accurate  estimate  of  x(t), 
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i.e.,  trace  P(t)  is  small,  the  weighting  factor  Q(t)P~*(t)  will  tend  to  be  "large*'. 
This,  of  course,  emphasizes  the  importance  of  this  second  correction  term  as  it 
should.  Indeed,  as  P(t)  0,  we  have  x(t)  *+  x(t),  and  we  would  expect  the  correc¬ 
tion  term  to  dominate  the  filter's  behavior  in  an  effort  to  force  x(t|t  +  T)  into 
correspondence  with  x(t). 

Let  us  now  relate  the  problem  that  we  have  solved  here  to  the  classical  Wiener 
smoothing  problem^**^ .  In  the  classical  formulation,  we  consider  the  block  dia¬ 
gram  shown  below  in  Fig.  2.  The  message  and  measurement  models  are  described  by 
Eqs.  (1)  and  (2),  respectively,  along  with  the  corresponding  statistical  information 
which  was  given  in  Section  1.0.  The  ideal  smoothing  filter  is  characterized 
by  its  n  x  n  system  Impulse  response  matrix  which  is 


i 
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At(t,T)  - 


t  *  T 

t  <  T 


j 

* 


) 

where  I  is  the  n  x  n  identity  matrix,  6  is  the  Dirac  delta  function,  and  T  >  0.  By 

the  notation  A^(t,T),  we  mean  the  response  or  state  of  the  filter  at  time  t  for  a  vector 
unit  impulse  (Dirac  delta  function)  input  at  time  T.  From  the  familiar  convolution 

integral, 

. 

i(t)  -  J  A^t.T)  x(T)  d  T  t  a  to  +  T 

fco 


J 

J 


I 


*  J  X(T)  6(t  -  T  -  T)  d  T 

fcO 

-  x(t  -  T) 


i 

V 
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i.e.,  the  ideal  state  i(t)  is  the  actual  state  delayed  by  T  units  of  time. 

Although  the  ideal  smoothing  filter  is  physically  realizable  by  a  pure  delay 
T,  the  operation  that  it  is  to  perform  cannot  be  implemented  since  we  do  not  have 
physical  access  to  x(t).  Hence,  it  is  the  "ideal"  against  which  we  compare  the 
performance  of  the  actual  smoothing  filter  which  operates  on  the  measurements.  We 
require  that  this  latter  filter  be  linear  and  that  its  n  x  m  impulse  response  matrix 
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A(t,T)  =  0  for  all  t  <  t  . 

Since  the  actual  smoothing  filter  requires  x(t^| t^  +  T)  as  its  initial  con¬ 
dition,  it  cannot  begin  its  operation  until  measurements  over  the  interval 
[t^jt^  +  T]have  been  processed  to  obtain  this  initial  condition.  Hence,  in 
smoothing,  we  require  that  t  ^  t^  +  T. 

Again,  utilizing  the  convolution  integral,  we  can  write 


x(t 


-  T|t)  =  J  A(t  -  T,T)  z(T)  d  T 


where  t  ^  tQ  +  T,  or  equivalently,  by  a  simple  change  of  variable,  we  also  have 


5t+T 

A(t,T)  z(T)  d  T 
t 


(44) 


where  now  t  ^  t  . 

o 

The  smoothing  error  is  defined  as 


x(t  -  T| t)  =  x(t  -  T)  -  x(t  -  Tf  t) 


for  t  ^  t^  +  T,  or,  equivalently. 


x(t|t  +  T)  =  x(t)  -  x(t|t  +  T) 


where  t  ^  t  . 

o 

We  then  say  that  a  filter  of  the  form  described  in  Eq.  (44)  that  minimizes  the 
mean-square  smoothing  error 

S  =  E[  x'(t|t  +  T)  x( 1 1 1  +  T)] 

for  all  t  ^  t  is  an  optimal  fixed-lag  smoothing  filter. 

In  the  classical  formulation,  we  see  that  the  problem  is  that  of  specifying 
the  filter  impulse  response  matrix  A(t,T).  The  results  obtained  in  this  paper  ex¬ 
press  the  differential  equation  for  the  filter  rather  than  its  impulse  response 
matrix.  Hov;ever,  the  latter  could  be  determined  from  Eqs.  (27)  and  (28)  if  de¬ 
sired. 

In  conclusion,  we  remark  that  the  question  of  the  stability  of  the  optimal 
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fixed-lag  smoothing  filter,  Eq.  (27),  its  gain  matrix,  Eq.  (37),  and  the  error 
covariance  matrix,  Eq.  (43),  remains  as  a  problem  area  for  future  study.  We  con¬ 
jecture  that  if  the  message  and  measurement  process  of  Eqs.  (1)  and  (2)  is  uniform¬ 
ly  completely  controllable  and  uniformly  completely  observable,  then  the  smoothing 
filter  is  uniformly  asmyptotically  stable  and  the  gain  and  covariance  equations  have 
equilibrium  solutions  for  t  *  -  00 . 
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